library(foreign)
library(cem)
library(MatchIt)
library(Zelig)
library(car)
library(Amelia)

setwd("~/Desktop/Replication")

stata <- read.dta("datafromraw.dta")

treatment <- matrix(data = stata$lrexclpop)
stata <- cbind(stata, treatment)
stata$treatment <- recode(stata$lrexclpop,"0 = 0;0.03:hi= 1")

######### Wimmer Model 1 and 2 Variables #############

vars1 <- c("newonset","newethonset","onsetstatus","actoraim", "lrexclpop", "egipgrps", "pimppast", "gdpcapl", "lpopl", "ongoingwarl", "year", "npeaceyears", "nspline1", "nspline2", "nspline3","cowcode","treatment")

stata1 <- stata[vars1]

vars2 <- c("newonset","newethonset","onsetstatus","actoraim", "lrexclpop","egipgrps", "pimppast", "ethfrac","gdpcapl", "lpopl", "ongoingwarl", "lmtnest", "regchg3", "anocl", "oilpcl", "year", "npeaceyears", "nspline1", "nspline2", "nspline3","treatment")

stata2 <- stata[vars2]

###############
#Missing Data
###############

set.seed(777)

bds <- matrix(c(13,0,150), nrow = 1, ncol = 3)

vari <-c("year","gdpcap","popavg","lpop","oilpc","anoc","democl","mtnest", "cowcode","newonset","newethonset","treatment","gdpcapl","lpopl","oilpcl","anocl","regchg3","ethfrac","eeurop","western","lamerica","ssafrica","asia","lmtnest","pimppast","egipgrps","lrexclpop","ongoingwarl","npeaceyears","nspline1","nspline2","nspline3","exclgrps","exclpop","ttlpop","discpop","pwrlpop","maxexclpop","intensity")

statai <- stata[vari]

imputed.sets <- amelia(statai, m = 5, ts = "year", cs = "cowcode", log = c("gdpcapl","egipgrps"))

imputed <- imputed.sets$imputations

################
# Recreate Wimmer Model 1 with multiply imputed data
################

output <- zelig(newonset ~ lrexclpop + pimppast + egipgrps + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3, model = "logit", data = imputed)

summary(output)

######
# Model showing insignificant interaction between lrexclpop and gdpcapl
######

output2 <- zelig(newonset ~ lrexclpop + pimppast + egipgrps + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3 + lrexclpop*gdpcapl, model = "logit", data = imputed)

summary(output2)


################       MODEL DIAGNOSTICS     ###############

################
# Matching on Treatment
################

# Find imbalance

imbalance(group = stata$treatment, data = stata[vars1], drop = c("newonset","newethonset","onsetstatus","actoraim","lrexclpop","nspline1","nspline2","nspline3","cowcode"))

### 0.0% local common support, L1 = 1 ###

########################################
# Nearest Neighbor Matching         ####
########################################

statai <- na.omit(statai)

nearest.match1 <- matchit(formula = treatment ~ egipgrps + pimppast + gdpcapl + lpopl + ongoingwarl + year + npeaceyears + nspline1 + nspline2 + nspline3, dataframe = imputed, data = statai, method = "nearest", distance = "logit")

nearest.data1 <- match.data(nearest.match1)

imbalance(group = nearest.data1$treatment, data = nearest.data1[vars1], drop = c("newonset","newethonset","onsetstatus","actoraim","lrexclpop","nspline1","nspline2","nspline3","cowcode"))

####################
### We cut a lot of "treatment" cases, but we still have L1 = 1!####################

output.nn1 <- zelig(newonset ~ lrexclpop + pimppast + egipgrps + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3, model = "logit", data = nearest.data1)

summary(output.nn1)

################
### lrexclpop remains significant, but gdpcapl remains much more significant ###
################


################ CEM ###################
########################################

gdpcut <- c(0,3.25,7.67,111)
popcut <- c(0,8.21,10,15)
npeaceyearscut <- c(0,1,3,11,100)
groupscut <- c(0,1,2,4,15)
yearcut <- c(1946,1960,1970,1980,1995,2005)

mat1 <- cem(treatment = "treatment", data = stata1, drop = c("newonset","newethonset","onsetstatus","actoraim", "lrexclpop","cowcode","nspline1", "nspline2","nspline3"), cutpoints = list(gdpcapl = gdpcut, npeaceyears = npeaceyearscut,  lpopl = popcut, egipgrps = groupscut,year = yearcut))

est1 <- att(mat1, newonset ~ lrexclpop + egipgrps + groups + gdpcapl + lpopl + year + ongoingwarl + year + nspline1 + nspline2 + nspline3, data = stata, model="logit")

summary(est1)

mat1

#### We lose many cases with only a marginal improvement in balance (L1 = 0.998!). None of our vars of central interest are significant, although note that gdpcapl survives very, very robustly). ####

############# WIMMER ET AL MODELS: FIRST DIFFERENCE ESTIMATIONS ############

x.exclpop.1 <- setx(output,lrexclpop=0)
x.exclpop.3 <- setx(output, lrexclpop= (3.3))
s.out.excl <- sim(output, x = x.exclpop.1, x1 = x.exclpop.3)

summary(s.out.excl)

x.gdp.1 <- setx(output, gdpcapl = 1.29)
x.gdp.3 <- setx(output, gdpcapl = 7.55)
s.out.gdp <- sim(output, x = x.gdp.3, x1 = x.gdp.1)

summary(s.out.gdp)

x.exclpop.lo <- setx(output,lrexclpop=0)
x.exclpop.hi <- setx(output, lrexclpop= 4.6)
s.out.exclmax <- sim(output, x = x.exclpop.lo, x1 = x.exclpop.hi)

summary(s.out.exclmax)

x.gdp.lo <- setx(output, gdpcapl = 7.6)
x.gdp.hi <- setx(output, gdpcapl = 21.6)
s.out.gdpmax <- sim(output, x = x.gdp.hi, x1 = x.gdp.lo)

summary(s.out.gdpmax)

x.gdp.lo2 <- setx(output, gdpcapl = 21.6)
x.gdp.hi2 <- setx(output, gdpcapl = 110)
s.out.gdpmax2 <- sim(output, x = x.gdp.hi2, x1 = x.gdp.lo2)

summary(s.out.gdpmax2)

x.extrapolate1 <- setx(output, gdpcapl = 21.6, lrexclpop = 4.6)
x.extrapolate2 <- setx(output, gdpcapl = 7.6, lrexclpop = 4.6)
s.out.gdpmax2 <- sim(output, x = x.extrapolate1, x1 = x.extrapolate2)

summary(s.out.gdpmax2)

####################

x.egipgrps.1 <- setx(output,egipgrps=1)
x.egipgrps.3 <- setx(output, egipgrps=2)
s.out.egipgrps <- sim(output, x = x.egipgrps.1, x1 = x.egipgrps.3)

summary(s.out.egipgrps)


x.egipgrps.lo <- setx(output,egipgrps=0)
x.egipgrps.hi <- setx(output, egipgrps= 14)
s.out.egipgrpsmax <- sim(output, x = x.egipgrps.lo, x1 = x.egipgrps.hi)

summary(s.out.egipgrpsmax)

################ REGIONAL INTERACTION MODELS ########################

### Fixed Effects ###

output.fe <- zelig(newethonset ~ lrexclpop + egipgrps + pimppast + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3 + factor(cowcode), model = "relogit", data = imputed)

summary(output.fe)

### lrexclpop interactions ###

output.re <- zelig(newethonset ~ lrexclpop*western + lrexclpop*eeurop + lrexclpop*lamerica + lrexclpop*ssafrica + lrexclpop*asia + lrexclpop + egipgrps + pimppast + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3 + lamerica + ssafrica + western + asia + eeurop, model = "relogit", data = imputed)

summary(output.re)

### Wimmer et al model with relogit and newethonset ###

output.old <- zelig(newethonset ~ lrexclpop + egipgrps + pimppast + gdpcapl + lpopl + ongoingwarl + npeaceyears + year + nspline1 + nspline2 + nspline3, model = "relogit", data = imputed)

summary(output.old)

####################
# First differences for interaction models, East Europe v Mid East
####################

### 3rd to Median exclusion

x.eeurop.2 <- setx(output.re,lrexclpop = 1.96, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.eeurop.3 <- setx(output.re, lrexclpop = 3.30,eeurop = 1,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out32 <- sim(output.re, x = x.eeurop.2, x1 = x.eeurop.3)

summary(s.out32)

x.nafrme.2 <- setx(output.re,lrexclpop = 1.96, eeurop = 0, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrme.3 <- setx(output.re, lrexclpop = 3.30,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.na32 <- sim(output.re, x = x.nafrme.2, x1 = x.nafrme.3)

summary(s.out.na32)

### Median to no exclusion

x.eeurop.lo <- setx(output.re,lrexclpop = 0, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.eeurop.hi <- setx(output.re, lrexclpop = 1.96,eeurop = 1,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out <- sim(output.re, x = x.eeurop.lo, x1 = x.eeurop.hi)

summary(s.out)

x.nafrme.lo <- setx(output.re,lrexclpop = 0, eeurop = 0, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrme.hi <- setx(output.re, lrexclpop = 1.96,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.na <- sim(output.re, x = x.nafrme.lo, x1 = x.nafrme.hi)

summary(s.out.na)

### Max to 3rd Quartile Exclusion

x.eeurop.max <- setx(output.re,lrexclpop = 4.6, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out4 <- sim(output.re, x = x.eeurop.3, x1 = x.eeurop.max)

summary(s.out4)

x.nafrme.max <- setx(output.re,lrexclpop = 4.6, eeurop = 0, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.na4 <- sim(output.re, x = x.nafrme.3, x1 = x.nafrme.max)

summary(s.out.na4)

###

x.eeurop1 <- setx(output.re,lrexclpop = 0, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrme1 <- setx(output.re, lrexclpop = 0,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.naee <- sim(output.re, x = x.eeurop1, x1 = x.nafrme1)

summary(s.out.naee)

x.eeurop <- setx(output.re,lrexclpop = 1.96, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrme <- setx(output.re, lrexclpop = 1.96,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.eena <- sim(output.re, x = x.eeurop, x1 = x.nafrme)

summary(s.out.eena)

x.eeurop3 <- setx(output.re,lrexclpop = 3.3, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrme3 <- setx(output.re, lrexclpop = 3.3,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.eena3 <- sim(output.re, x = x.eeurop3, x1 = x.nafrme3)

summary(s.out.eena3)

x.eeuropmax <- setx(output.re,lrexclpop = 4.6, eeurop = 1, asia = 0, western = 0, lamerica = 0, ssafrica = 0)
x.nafrmemax <- setx(output.re, lrexclpop = 4.6,eeurop = 0,asia = 0, western = 0, lamerica = 0, ssafrica = 0)
s.out.eenamax <- sim(output.re, x = x.eeuropmax, x1 = x.nafrmemax)

summary(s.out.eenamax)